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An efficient all-electron G^W" method and a quasiparticle selfconsistent GW (QSGW) method for molecules 
are proposed in the molecular orbital space with the full random phase approximation. The convergence with 
basis set is examined. As an application, the ionization energy (/) and electron affinity (A) of a series of 
conjugated molecules (up to 32 atoms) are calculated and compared to experiment. The QSGW result improves 
the G'^W" result and both of them are in significantly better agreement with experimental data than those from 
Hartree-Fock (HF) and hybrid density functional calculations, especially for A. The nearly correct energy gap 
and suppressed self-interaction error by the HF exchange make our method a good candidate for investigating 
electronic and transport properties of molecular systems. 



Understanding and manipulating the electron transport 
through molecules is the basis of molecular electronics which 
provides a promising way for achieving the ultimate miniatur- 
ization of electronic components |1J. To be enable to calcu- 
late accurately the transport properties of molecular devices 
is therefore desirable for the interests of physics and appli- 
cations. In this regard, a standard method which has been 
extensively used is the density functional theory (DFT) |l2l] 
combined with the nonequilibrium Green function (NEGF) 
method |3]. Despite the success of DFT in describing ground- 
state electronic properties, its application to quantum trans- 
port faces much more challenges |4-6]. Intrinsic issues with 
DFT, such as self-interaction error (SIE) and underestimated 
energy gap, will lead to incorrect molecule-lead charge trans- 
fer, incorrect position of the chemical potential in the gap be- 
tween the highest occupied molecular orbital (HOMO) and 
the lowest unoccupied molecular orbital (LUMO) as well 
as the broadening of the HOMO and LUMO, especially for 
weakly coupled systems 1^. Consequently, the molecular 
conductance will be overestimated, even by orders of mag- 
nitude in some cases, as has been shown by some previous 
work in literature Si ■ 

A good electronic structure method for a NEGF-based 
transport calculation should satisfy several conditions: (i) 
nearly correct HOMO and LUMO energies and the gap, (ii) 
suppressed SIE, (iii) conserving the number of particle (CNP), 
and (iv) implemented with localized basis functions. Physi- 
cally, even an exact DFT calculation with a local exchange- 
coiTelation potential (Vxc) cannot give a coiTect HOMO- 
LUMO gap due to the lack of the discontinuity. Improve- 
ment with a nonlocal Vxc is still under development and fac- 
ing challenges |6-8]. Different from DFT, Hartree-Fock (HF) 
method is SIE free and can usually give a reasonable HOMO 
energy (Koopermann's theorem) but lacks electron coiTelation 
and cannot predict a reasonable LUMO energy. Quasi-particle 
calculation in the GW approximation |9-11] may provide a 
solution but its non-selfconsistent version (G'^W'^) has a crit- 
ical issue with CNP which has been shown to be important 
for transport calculations Ml2l 11311 . It also faces challenging 
computational effort for large molecular systems. Addition- 
ally, in spite of its very successful applications for bulk semi- 



conductors] 14- 1611. it accuracy for molecules is still an open 
problem |fTlUl7m9ll considering the stronger electronic re- 
laxation effect ilTL Furthermore. as shown in a recent work 
1 11], the core electrons play an important role in determining 
accurately the excitation energies of molecules. Therefore, an 
efficient ab initio all-electron GW approach being able to deal 
with large molecules is desirable. 

As an effort towards a more accurate quantum transport cal- 
culation, in this paper, we first propose an efficient all-electron 
non-selfconsistent G"W" method based on the full random 
phase approximation (RPA) using Gaussian basis functions 
and a HF input. To achieve a high computational efficiency, it 
is implemented in the molecular orbital (MO) space with tech- 
niques for reducing the error coming from the incompleteness 
of the basis set. The correlation self-energy is determined 
first on the imaginary energy axis and then on the real en- 
ergy axis using the analytical continuation approach which 
was proposed originally for a space-time approach ll20i I21I1 . 
The convergence with respect to the size of basis set is exam- 
ined by calculating He and Be atoms using a series of basis 
sets. 

To address the issues with the G^W*^ method, we fur- 
ther implement a q uasiparticle selfconsistent GW (QSGW) 
method I22II . in which a selfconsistency is performed 

between the G"W° and a quasi-DFT calculation with a non- 
local Vxc constructed from the G^W'^ selfenergy. In this way, 
the converged Vxc will largely suppress the SIE since the real 
(nonlocal) exchange is used and is also CNP compliant. Fur- 
thermore, as is shown in this work, the quasi-DFT calculation 
with this Vxc can give a nearly correct HOMO and LUMO en- 
ergies, even better than the G^W^ result, implying that it will 
satisfy the four conditions for a good NEGF-based quantum 
transport calculation. 

As an application of our methods, we investigate the ion- 
ization energy (/) and electron affinity (A) of a series of con- 
jugated molecules. As is well known, conjugated molecules 
play the most important roles in molecular electronics be- 
cause of their small HOMO-LUMO gap, good molecular 
conductance, and tunable electronic states. The systems 
studied include acetylene (C2H2), ethylene (C2H4), allene 
(C3H4), diacetylene (C4H2), benzene (CeHg), phenylacety- 
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lene (CsHg), naphthalene (CioHg), biphenyl (C12H10), an- 
thracene (C14H10), and perylene (C20H12), for which rehable 
experimental data are available 02311 . Comparison to the ex- 
perimental data shows that our G^W*^ and QSGW results im- 
prove the HF ones significantly, especially for A. The QSGW 
improves the G^W^ further and shows a very good agree- 
ment between theory and experiment. This indicates that the 
all-electron GW calculation can describe very well molecu- 
lar electronic structures. The computational efficiency of our 
methods makes it possible to easily deal with systems con- 
sisting of several tens of atoms from the first principles. The 
satisfaction of the four conditions for quantum transport im- 
plies the QSGW-based quasi-DFT approach a good candidate 
for investigating electronic and transport properties of molec- 
ular devices. 

The GW selfenergy is divided into a bare-exchange part and 
a correlation part: S = + S^. On the imaginary energy 
axis the correlation part is 

1 r°° 

= -— / Gir,r',tiLu + cj'))W%r,r';iu;')dw\(l) 
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where G is the Green function and = W — v is the 
screened Coulomb potential minus the bare Coulomb poten- 
tial. In the MO space (the MOs are denoted by m, n, ...) the 
Green function simply reads, 



(2) 



lUJ - Cn 

with e„ denoting the nth eigenvalue. 

W''{r,r';iu}) = [5{r ~ r ) - P{r,r';iuj)v]~^ v - v,(3) 
and we define 

W%r,r';iuj) = [5(r - r') - P(r, r'; - 1, (4) 

where P is the polarization function which can be determined 
by RPA: 
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where / is occupation number and i/j's are the MOs. If the 
MO space is complete, P can be expressed as 



^ — ' £; — ei- — VjJ 

kl 

in terms of the three-center overlap integrals of the MOs, 



(6) 
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However, in practice, the incompleteness of the MO space will 
cause error in the calculation of P, and even worse, causes 
larger error in the product of Pv. To suppress this error, we 



avoid the individual calculations of P and v but calculate the 
product as a whole in terms of the electron integrals, 

[^(-)-]™. = E , Q^c^', (8) 
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where. 



(9) 



are the three-center Coulomb integrals of the MOs. The two 
kinds of three-center MO integrals are determined in terms 
of the corresponding three-center integrals of the atomic or- 
bitals (AO) presented by the standard contracted Gaussian ba- 
sis functions, which can then be calculated analytically. After 
combining Eqs. (1), (2), (4), and (8), we get 
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which is simply a product of the three three-dimensional ar- 
rays. The integral Wjki can be calculated by Gaussian quadra- 
ture. Because W varies extreme smoothly on the imaginary 
energy axis, only very few energy grids (about 100 points 
from to 00) are needed to well converge the elements of 
W. 

After having obtained S'^ on the imaginary energy axis, we 
determine its values on the real energy axis using the analyt- 
ical continuation approach ||20[|2lll . and least-square fit each 
element to the multipole form using a BEGS optimizaton tech- 
nique, oo + ui-b- ^''■^ complex parameters and bi. 
For bulk semiconductors it was found that n 2 is usually 



good enough 11201 12111 . For molecules, however, we find that 
n = 6 is generally required to give an extremely stable and 
accurate fitting. 

To address the issues with the G'^W'^ method, we imple- 
ment a QSGW approach, in which the G'^W'^ selfenergy is 
used to construct a nonlocal exchange-correlation potential 
Vxc for a quasi-DFT calculation. In practice, we follow the 
scheme proposed in which has been shown 

to be very successful for bulk semiconductors, but now imple- 
mented in the MO space for molecules. 



[V^cUn = Re + S^(en)]„„, /2 + 



(11) 



with being calculated in terms of the HF exchange oper- 
ator The quasi-DFT calculation will generate new eigenval- 
ues and MOs which can be used for the next G°W° calcula- 
tion. This procedure is going on until the maximum change 
in the Vxc elements is smaller than the convergence criteria. 
Previously, the QSGW approach implemented in plane-wave 
pseudopotential formalism was applied to bulk semiconduc- 
tors and improved significantly the G°W'^ band gaps 1 1412211 . 
For molecules, its applicability and quality is still unknown. 



3 



26.0 
25.8 
25.6 
25.4 
25.2 
25.0 
24.8 
24.6 
24.4 
24.2 





' 








1 1 1 

He: HOMO 
— HF 
-©-qOw" 

Expt 














































2 


5 


3 1 


5 35 55 70 126 182 



Number of basis function 





































3 C 


Be: HOMO 

-O-HF 

-m-G°w° 

Expt 


1 1 

> ( 




J C 










1 


5 1 


5 19 35 70 105 126 



Number of basis function 

FIG. 1 : Dependence of the -HOMO energy of He atom (upper panel) 
and Be atom (lower panel) on the size of basis set used in the HF and 
G^W" calculations. The experimental result is shown by a horizontal 
dashed line. 



Like other correlated electronic structure methods (for in- 
stance, the second-order M0ller-Plesset theory (MP2) and 
coupled-cluster method) the G^W^ method is more sensitive 
to the size of basis set than DFT and HF method. So far, care- 
ful examinations for its convergence behavior are still lack- 
ing in literature llllll . probably leading to the scattering re- 
sults from different calculations and biased conclusions for 
different systems. To examine the basis-set convergence of 
our method, we first calculate the -HOMO energy (i.e., /) of 
He and Be atoms, adopting different basis sets l24ll : 6-31G*, 
6-31 IG**, cc-pvDZ, cc-pvTZ, cc-pvQZ, aug-cc-pvQZ, cc- 
pv5Z, cc-pv6Z, and aug-cc-pv6Z, and plot the results in Fig. 
1. 

For He atom, the 6-31G*, cc-pvDZ, and 6-311G** basis 
sets have only 2, 5, and 6 basis functions, respectively. For 
such small basis sets its G^W'^ -HOMO energy fluctuates re- 
markably. However, once the basis set is larger than 15 basis 
functions the results of both He and Be atoms tend to be well 
converged and the fluctuation is around O.leV. The important 
thing is that for both atoms the G'^W'^ result improves signif- 
icantly the HF result: For He the -HOMO energy is corrected 
downwards while for Be it is coiTected upwards, becoming 
much closer to the experimental data. 

For a molecule, its I or A includes two contributions, i.e., 
the vertical / or A and the corresponding charge-induced 
structural relaxation energy (AE^ or AE~): I = /vert. - 
AE'^ and A = Avert. + AE^. AE+'~ can be determined by 
performing a structural relaxation calculation with the charge. 
For the ten conjugated molecules studied in this work we first 



optimize their structures on the DFT/B3LYP/cc-pvTZ level, 
and then calculate their Ai?+ and AE^ by allowing further 
relaxation with the charge. From the results listed in Table 
1 one can see that AE^ is comparable with AE^ and their 
values depend strongly on the molecular structure. 

The results of / and A of the ten conjugated molecules 
given by the different methods are listed in Table 1 together 
with the available experimental data cited from Ref.|23]. For 
A only the four largest molecules are studied, which have a 
positive or nearly zero experimental value. The results are 
also ploted in Fig. 2 for a better vision and comparison. For 
molecular systems the hybrid DFT/B3LPY has been proven to 
be a significant improvement to DFT/LDA and DFT/GGA be- 
cause the HF exchange is partly included, suppressing partly 
the SIB with LDA and GGA 0. However, as shown in Table 
1 and Fig. 2, its results of / and A of the conjugated molecules 
are still far away from being correct. In spite of lacking elec- 
tron correlation, HF method gives much better result of / due 
to the Koopermann's theorem, but its result of A is very bad, 
being comparable to the B3LYP result but with an opposite 
trend (see Fig. 2). The total-energy difference calculation us- 
ing HF (AHF) improves HF for A considerably but AHF re- 
sult of / becomes much worse because of the inaccurate total 
energy. On the other hand, the total-energy difference calcula- 
tion using B3LYP (AB3LPY) gives very good results of both 
/ and A because of its accurate total energy. 

The G'^W" result improves the HF result considerably for / 
while significantly for A, showing that the correlation is more 
important for the excited states. The QSGW result improves 
the G°W° result further, being in very good agreement with 
the available experimental data [23] and the AB3LPY result. 
This finding is consistent with previous calculations fl^, [3. 



2211 which showed a significant improvement from QSGW to 
G^W" for band gaps in bulk semiconductors. We note that the 
quasi-DFT calculation with the nonlocal Vxc from the QSGW 
method can be a good candidate for NEGF-based quantum 
transport calculations considering the facts that (i) it can give 
nearly correct HOMO and LUMO energies of molecules, (ii) 
it is CNP compliant, (iii) SIE is largely avoided because of 
the use of the HF exchange, and (iv) it is implemented with 
the local Gaussian basis functions. 

Finally, our work shows that the all-electron ab initio G'^ W'^ 
and QSGW methods can describe very well the electronic 
structures of conjugated molecules. This finding is consis- 
tent with a previous ab initio G°W" calculation for small 
molecules ll ill , which also found that the inclusion of core- 
valence exchange is very important in obtaining accurate exci- 
tation energies, but is somewhat different from a recent 7r-only 
model GW calculation IlITll which shows that the accuracy of 
GW method is largely affected by the electron-removal and 
electron-addition induced relaxation. This difference may be 
due to the role played by the a states which have much broader 
energy range, and the resulting tt-ct relaxation. 

This work was supported by Shanghai Pujiang Pro- 
gram under Grant 10PJI4I0000 and the MOST 973 Project 
2011CB922200. 
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TABLE I: Calculated (cc-pVTZ level) and experimental results of I and A (in eV) of ten conjugated molecules, as well as their charge-induced 
relaxation energies AE^ and AE^ . The experimental data are cited from Ref. 12311 . 
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FIG. 2: Comparison between the calculated results and experimental data of the / and A of the ten conjugated molecules in Table 1. 
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